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THE  GENERALIZED  SUNDMAN 
TRANSFORMATION  FOR  PROPAGATION  OF 
HIGH-ECCENTRICITY  ELLIPTICAL  ORBITS 

Matthew  Berry*  Liam  Healy^ 


Abstract 

A  generalized  Sundman  transformation  dt  =  cr^ds  for  exponent  n  >  1  may  be  used  to 
accelerate  the  numerical  computation  of  high-eccentricity  orbits,  by  transforming  time  t  to  a 
new  independent  variable  s.  Once  transformed,  the  integration  in  uniform  steps  of  s  effectively 
gives  analytic  step  variation  in  t  with  larger  time  steps  at  apogee  than  at  perigee,  making  errors 
at  each  point  roughly  comparable.  In  this  paper,  we  develop  techniques  for  assessing  accuracy 
of  s-integration  in  the  presence  of  perturbations,  and  analyze  the  effectiveness  of  regularizing 
the  transformed  equations.  A  computational  speed  comparison  is  provided. 


INTRODUCTION 

Naval  Space  Command  (NSC)  has  adopted  a  special  perturbations  catalog  system,  that  is,  orbit 
propagation  based  on  a  numerical  integrator.  It  is  more  accurate  than  a  general  perturbations  system 
which  is  based  on  analytic  theory.  The  integrator  forms  a  core  routine  in  the  parallel-processing  space 
surveillance  computation  suite  of  programs  known  as  SpecialK  (Ref.  1).  It  is  used  for  ephemeris 
generation  and  differential  correction  by  both  the  automated  parallel  program  and  the  manual 
maintenance  program.  This  suite  is  used  to  maintain  approximately  1300  satellites  of  greatest 
interest  for  space  surveillance  work;  shortly,  that  number  of  satellites  will  grow  to  approximately 
3300  and  then  the  current  catalog  of  10,000  objects  will  be  processed  this  way.  Within  the  next 
ten  years  the  NSC  Service  Life  Extension  Program  will  produce  observations  on  50,000  -  100,000 
objects;  these  orbits  will  be  maintained  using  special  perturbations. 

Because  of  the  large  number  of  satellites,  attention  must  be  paid  to  the  total  time  of  computation. 
For  very  eccentric  orbits,  satisfactory  accuracy  can  be  maintained  with  a  larger  step  size  near  apogee 
than  is  used  at  perigee.  One  way  to  use  a  larger  step  at  apogee  is  to  make  use  of  the  known  two-body 
orbit  to  transform  the  independent  variable  from  time  t  to  another  variable  s,  dt  =  cr^ds.  This 
transformation,  a  generalized  Sundman  transformation  (Ref.  2),  allows  one  to  use  any  numerical 
integration  method  in  the  independent  variable  s,  which  is  thus  called  s-integration.  Conventional 
time  integration,  by  contrast,  is  called  t-integration. 

This  paper  presents  a  description  of  s-integration,  including  implementation,  as  well  as  a  dis¬ 
cussion  and  analysis  of  those  orbits  for  which  s-integration  is  preferable  to  t-integration.  Since  the 
goal  of  using  s-integration  is  to  save  computation  time,  relative  speed  comparisons  of  s-  versus  t- 
integration  are  presented  here.  For  orbits  of  sufficient  eccentricity,  the  comparisons  show  a  dramatic 
speedup  over  fixed-step  ^integration,  thus  validating  the  original  motivation  for  using  s-integration. 

*  Graduate  Assistant,  Department  of  Aerospace  and  Ocean  Engineering,  Virginia  Tech,  Blacksburg,  Virginia  24061, 
and  Naval  Research  Laboratory,  Code  8233,  Washington,  DC  20375-5355,  E-mail:  maberry2@vt.edu. 

^Research  Physicist,  Naval  Research  Laboratory,  Code  8233,  Washington,  DC  20375-5355,  and  Lecturer,  Depart¬ 
ment  of  Aerospace  Engineering,  University  of  Maryland,  College  Park,  MD  20742.  E-mail:  Liam.Healy@nrl.navy.mil. 
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We  conclude  with  eccentricity  values  above  which  s-integration  is  advisable  because  there  is  a  speed 
advantage  without  loss  of  accuracy. 

GENERALIZED  SUNDMAN  TRANSFORMATIONS 

Transformation  equation 

Sundman  (Ref.  3)  and  Levi-Civita  (Ref.  4),  in  attempting  to  solve  the  restricted  problem  of  three 
bodies,  introduced  the  transformation  of  the  independent  variable 

dt  =  cr  ds,  (1) 

with  c  constant  for  the  two-body  orbit,  because  this  transformation  regularizes,  and  in  fact  linearizes, 
the  equations  of  motion.  Later  investigators  raised  r  to  different  powers  in  the  transformation, 

dt  =  cr'^  ds,  (2) 

known  as  the  generalized  Sundman  transformation  (Szebehely  and  Bond  (Ref.  5)  generalized  even 
more,  allowing  an  arbitrary  function  of  r). 

The  new  independent  variable  s  is  better  understood  in  terms  of  an  orbit  angle.  We  give  this 
term  to  any  angle  9  considered  a  function  of  true  anomaly  9{u)  that  has  the  following  properties: 

•  At  perigee,  the  value  of  the  angle  is  the  same  as  true  anomaly:  9{u)  =  u  =  2'Km  for  any  integer 

m. 


•  At  apogee,  the  value  of  the  angle  is  the  same  as  true  anomaly:  9{u)  =  n  =  nm  for  any  odd 
integer  m. 

•  The  angle  increases  monotonically  with  true  anomaly,  9{u2)  >  9{ui)  if  U2>  oi. 

•  There  is  symmetry  about  the  major  axis:  9{i')  =  —9{—i>). 

Examples  include  the  mean,  eccentric,  and  true  anomalies.  Simply  applying  the  transformation  to 
s  does  not  assure  that  s  is  then  an  orbit  angle;  c  must  be  picked  so  that  the  appropriate  boundary 
conditions  are  satisfied.  The  c  necessary  for  the  case  n  =  3/2  is  computed  in  detail  in  the  next 
section. 

For  each  of  the  possible  values  of  n  >  1,  there  is  a  corresponding  angle  (Ref.  6,  p.  100;  Ref.  7; 
Ref.  8,  p.  19;  Ref.  9,  p.  484): 

•  n  =  1  or  dt  =  cr  ds.  The  angle  s  is  the  eccentric  anomaly  if  c  is  chosen  so  that  s  is  an  orbit 
angle,  c  =  ^Jaj [i.  This  case  is  Sandman’s  original  transformation,  or  the  Kustaanheimo-Stiefel 
transformation  (Ref.  8). 

•  n  =  3/2ordt  =  cr‘^l^ds.  We  shall  focus  on  this  transformation. 

•  n  =  2  or  dt  =  cr‘^ds.  The  angle  s  is  the  true  anomaly  if  c  is  chosen  so  that  s  is  an  orbit  angle, 

c=  [/ia(l-e2)]-i/2. 

The  constants  c  necessary  for  the  cases  of  n  =  1  and  n  =  2  can  be  derived  from  partial  derivatives 
of  the  orbit  angles.  The  partial  derivative  of  the  eccentric  anomaly  with  respect  to  the  mean  anomaly 
is 

—  =  -  (3) 

dM  r’  ^  ’ 

with  a  the  semimajor  axis.  Because  M  =  mt,  where  m  =  fi/a^  is  the  mean  motion,  one  may 

describe  the  differential  relationship  between  time  and  eccentric  anomaly  E, 

dt^  —dM  =  —dE  ^r.  SdE.  (4) 

m  am  y  fi 
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The  relation  between  time  and  the  true  anomaly  may  be  developed  as  follows.  The  partial  derivative 
of  true  anomaly  with  respect  to  the  mean  anomaly  is  easy  to  compute, 


diy  _ 

where  e  is  the  eccentricity,  M  is  the  mean  anomaly  and  u  is  the  true  anomaly  (Ref.  10).  This  relation 
gives  us  the  time  rate  of  change,  since  time  is  proportional  to  mean  anomaly.  Thus 

^2 

dt  =  — .  =di'.  (6) 

VMi  -  e") 

Ferrer  and  Sein-Echaluce  (Ref.  11)  showed  that  only  n  =  1  and  n  =  2  linearize  the  Kepler 
problem;  the  latter  only  with  regularization.  Palmer  et  al.  (Ref.  12)  studied  the  use  of  an  n  =  1 
transformation  for  integration  with  a  Bulirsch-Stoer  integrator  and  compared  results  using  positions 
obtained  from  a  GPS-equipped  satellite,  but  as  described  below,  n  =  3/2  is  more  commonly  used  for 
integration.  Since  perturbations  affect  the  orbit,  they  also  affect  these  transformations.  However, 
the  analysis  that  follows  assumes  a  Kepler  problem  in  order  to  understand  the  general  characteristics 
of  s-integration. 

There  are  two  ways  we  can  approach  a  numerical  integration  implementation  using  these  trans¬ 
formations.  We  can  determine  that  the  step  size  at  perigee  is  fixed  through  the  transformation;  in 
this  case,  the  total  number  of  steps  varies  depending  on  the  transformation.  Expressing  the  transfor¬ 
mations  as  orbit  angles,  Eigure  1  illustrates  this  approach  with  an  approximately  fixed  perigee  step 
of  Ah'  1.0  radian.  Alternatively,  one  can  fix  the  total  number  of  steps  on  an  orbit,  and  allow  the 
step  sizes  to  vary.  Eigure  2  illustrates  this  approach  with  sixteen  points  equally  distributed  in  the 
same  four  orbit  angles.  The  case  of  n  =  3/2  represented  as  an  orbit  angle,  called  the  intermediate 
anomaly,  is  discussed  in  the  next  section. 


n  =  3/2,  or  intermediate  anomaly 


Merson  (Ref.  13)  introduced  the  idea  of  using  the  value  n  =  3/2  in  the  generalized  Sundman  trans¬ 
formation,  with  an  intent  not  of  regularization  per  se,  but  of  maximizing  computational  efficiency; 
see  also  Ref.  14.  He  gave  an  analysis  showing  that  this  value  of  n  equally  distributes  the  integration 
error  around  a  full  orbit,  even  if  the  eccentricity  is  high.  One  may  conclude  that  this  method  is 
preferred  for  numerical  integration.  He  also  gave  accuracy  and  timing  results  for  n  =  3/2,  compared 
to  other  integrators.  Nacozy  (Ref.  7),  like  Merson,  expressed  s  in  terms  of  an  elliptic  integral  of  the 
true  anomaly,  and  dubbed  this  angle  the  intermediate  anomaly.  His  choice  of  constant  c  =  1/ 
made  s  dimensionless  but  does  not  result  in  an  orbit  angle;  we  show  how  to  make  s  an  orbit  angle 
in  this  section. 

Eor  any  exponent  n,  we  may  express  the  s  differential  element  in  terms  of  the  true  anomaly 
differential  element  using  (6), 


ds  = 


1 

- ,  =dy. 

c  y/ia(l  -  e^) 


(7) 


Eor  n  =  3/2,  this  relation  simplifies  to 


ds  = 


1  dv 
c^/JI  Vl  +  ecosi/’ 


(8) 


using  the  relationship  r  =  p/{l  ecosn).  Introducing  the  half  angle  9  =  i^/2,  the  relation  can  be 
rewritten  via  a  trigonometric  identity 


ds  = 


2 _ dO 

+  e)  sin^  9  ’ 


(9) 
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(a)  Equal  mean  anomaly  (n=0)  with  58  steps. 


(b)  Equal  eccentric  anomaly  (n=l)  with  16 
steps. 


(c)  Equal  intermediate  anomaly  (n=3/2)  with  (d)  Equal  true  anomaly  (n=2)  with  6  steps. 

10  steps. 

Figure  1:  Points  separated  by  equal  values  of  various  orbit  angles,  with  the  corresponding  exponent 
n.  Each  orbit  has  approximately  the  same  step  size  Aiy  =  1  radian  at  perigee. 
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(a)  Equal  mean  anomaly  (n=0);  step  size  at  (b)  Equal  eccentric  anomaly  (n=l);  step  size  at 

perigee  is  Au  =  1.97  radians.  perigee  is  Au  =  0.97  radians. 


(c)  Equal  intermediate  anomaly  (n=3/2);  step  (d)  Equal  true  anomaly  (n=2);  step  size  at 

size  at  perigee  is  Au  =  0.60  radians.  perigee  is  Au  =  0.39  radians. 

Figure  2:  Points  separated  by  equal  values  of  various  orbit  angles,  with  the  corresponding  exponent 
n.  There  are  16  steps  in  each  orbit. 
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where  k‘^  =  2e/(l  +  e).  Integrating  (9),  the  second  fraction  is  the  incomplete  elliptic  integral  of  the 


first  kind  F{9,  k) 

S-  ^ 

+  ^2  / 

(10) 

In  order  to  find 
we  have  s(0)  = 
Substituting  the 

c,  we  must  match  the  boundary  conditions.  Evaluating  (10)  at  perigee  {y  =  0), 
0  because  F{0,k)  =  0.  Evaluating  (10)  at  apogee,  we  require  s  =  tt  at  y  =  tt. 
complete  elliptic  integral  of  the  first  kind  K{k)  (Ref.  15,  formula  17.3.2), 

TT-  Fl _ fFF-  Fl _ K[k), 

CA/p(l  +  e)  V2  )  m//i(l  +  e) 

(11) 

so 

2K{k) 

TT^/iJ.{l  +  e) 

(12) 

This  value  of  c  is  what  Merson  (Ref.  13,  p.  17)  called  an  “augmentation  factor.”  Outside  of  the 
half-plane  between  perigee  and  apogee,  0  <  <  tt,  the  elliptic  integral  should  be  evaluated  by 

making  use  of  the  elliptic  integral  magnitude  recursion  relation  (Ref.  15,  formula  17.4.3), 


Filn  ±  0,  m)  =  21K  ±  -F(0,  m). 


(13) 


With  ly  =  y  mod  27r,  the  change  in  true  anomaly  from  the  previous  perigee  passage,  and  N  =  [i^/27r], 
the  number  of  whole  orbits  completed,  s  may  be  expressed  as  an  orbit  angle. 


27riV  + 


if  0  <  1/  <  TT, 


27r(iV  + 1)  - (^,A:)  if  tt  <  P  <  27r. 


(14) 


This  definition  ensures  symmetry  about  the  major  axis. 

To  solve  for  the  true  anomaly  in  terms  of  the  intermediate  anomaly,  we  must  use  the  inverse  of 
the  elliptic  integral,  the  Jacobian  elliptic  function.  With 


we  may  write 


u 

sm  —  =  sn 

Zi 


and  in  the  general  case. 


V  = 


21^ N  +  2  arcsin  sn  (  — ^  s 
2'k{N  +  1)  —  2  arcsin  sn  g 


if  0  <  s  <  TT, 


if  TT  <  s  <  27r, 


(15) 

(16) 


(17) 


where  s  =  s  mod  27r  and  N  =  [s/27r].  These  definitions  complete  the  relationship  between  true 
anomaly  y  and  s  as  an  orbit  angle  when  n  =  3/2.  What  follows  in  the  rest  of  this  paper  does  not 
depend  on  s  being  an  orbit  angle  (that  is,  c  can  be  any  dimensionless  number  times 

Based  on  these  formulas,  it  is  straightforward  to  compute  the  ratio  of  the  number  of  steps  in  a 
complete  orbit  for  s-integration  and  t-integration  when  the  two  have  the  same  step  size  at  perigee 
(Fig.  1(a)  and  Fig.  1(c)).  First,  the  constant  step  size  in  a  given  orbit  angle  is  computed  by 
converting  the  value  in  true  anomaly  at  the  first  step  after  perigee  to  the  orbit  angle.  Because  all 
orbit  angles  are  zero  at  perigee,  this  first  step  is  the  step  size  at  perigee  in  the  orbit  angle.  For  the 
mean  anomaly  (t-integration)  the  conversion  uses  Kepler’s  equation;  for  the  intermediate  anomaly 
(n  =  3/2  s-integration),  we  use  (14).  Since  these  step  sizes  are  constant  throughout  the  orbit  for 
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the  respective  integration  methods,  the  ratio  of  the  number  of  points  is  the  inverse  of  the  ratio  of 
the  step  sizes.  This  ratio  is  plotted  as  a  function  of  eccentricity  in  Figure  3.  For  circular  orbits 
(e  =  0)  the  number  of  integration  steps  is  exactly  the  same;  as  eccentricity  increases,  t-integration 
has  proportionately  more  steps  than  s-integration.  One  might  be  tempted  to  conclude  that  it  is 
never  disadvantageous  to  use  s-integration  because  it  never  has  fewer  integration  steps  if  the  perigee 
step  size  is  chosen  the  same  as  t-integration.  However  this  conclusion  is  not  correct;  additional 
computational  costs  are  associated  with  s-integration  in  transforming  the  independent  variable  and 
integrating  the  transformation  equation  (2)  that  are  not  present  for  t-integration.  Furthermore, 
perturbations  affect  the  accuracy  of  s-integration  more  than  that  of  t-integration,  as  we  demonstrate 
in  later  sections. 


Ecxientricity 


Figure  3:  Ratio  of  the  integration  steps  per  orbit  for  t-integration  to  s-integration. 


IMPLEMENTATION  OF  THE  TRANSFORMATION 

SpecialK  has  an  eighth-order  Gauss-Jackson  integrator;  it  allows  the  user  to  choose  either  t-integration 
or  s-integration.  The  s-integration  is  the  Merson  /  Nacozy  form  discussed  in  the  previous  section 
with  n  =  3/2,  but  without  regularization,  and  with  c  =  l/y^,  so  s  is  not  an  orbit  angle.  The 
implementation  of  t-integration  is  described  in  detail  in  (Ref.  16).  Integration  in  s  uses  the  same 
code,  but  adds  the  necessary  transformation  from  t  to  s  space;  the  transformation  back  into  t  space 
requires  the  integration  of  an  additional,  seventh,  differential  equation  (2)  to  determine  the  time  t. 

The  software  first  selects  a  step  size  in  time.  For  s-integration  the  step  size  must  be  converted 
into  s-space,  and  the  conversion  is  performed  using  the  distance  from  earth  center  to  the  satellite  at 
perigee  Vp, 

As  =  ^/]2  Vp  ^  At.  (18) 

This  makes  the  step  at  perigee  the  same  in  s-integration  [Fig.  1(c)]  as  it  is  in  t-integration  [Fig.  1(a)]. 

After  the  step  size  is  determined  the  integrator  enters  its  initialization  phase,  which  is  necessary 
because  the  Gauss-Jackson  method  is  a  multi-step  integrator.  In  the  initialization  phase  the  8 
positions  and  velocities  surrounding  epoch  are  first  found  using  a  fifth-order  Taylor  expansion  of  the 
two-body  solution,  known  as  the  /  and  g  series  (Ref.  9)  (Eq.  4-68).  The  t-integration  initialization 
can  be  summarized: 


Initialization  for  t-integration 

1.  Use  /  and  g  series  to  calculate  8  positions  and  velocities  surrounding  epoch. 

2.  Evaluate  9  accelerations  from  the  positions  and  velocities,  including  epoch. 

3.  While  the  accelerations  have  not  converged: 

(a)  For  each  point  n  =  —4 ...  4,  n  7^  0: 

i.  Calculate  and  using  Gauss- Jackson  and  summed  Adams  mid-corrector  formu¬ 
las. 

ii.  Evaluate  the  acceleration  using  the  appropriate  force  model. 

(b)  Test  convergence  of  the  accelerations. 

A  more  complete  description  is  available  (Ref.  16,  pp.  18-19). 

Because  the  /  and  g  series  equations  depend  on  the  time  between  the  points,  and  for  s-integration 
the  points  must  be  equally  spaced  in  s,  a  conversion  must  be  made  from  s  to  time.  To  be  exact, 
the  time  should  be  found  by  integrating  (2);  however,  accuracy  is  not  as  critical  in  this  first  phase 
of  the  initialization  routine,  so  the  time  is  approximated  by  holding  the  epoch  distance  constant. 


After  the  positions  and  velocities  are  calculated  by  the  /  and  g  series,  the  force  model  is  evaluated 
to  find  the  accelerations.  These  accelerations  must  then  be  converted  to  s-space  by  changing  t 
derivatives  to  s  derivatives. 

The  derivatives  with  respect  to  s  are  computed  using  the  derivatives  with  respect  to  t  via  the 
relation 

4  =  c-v-"4.  (20) 

dt  ds 


Differentiation  with  respect  to  s  will  be  indicated  with  a  prime  (');  that  with  respect  to  t  by  a  dot 
(').  The  velocity  is  converted  by  application  of  (20), 


The  acceleration  can  then  also  be  transformed, 

r  =  -c-^c'r-^^r'  -  nc-\-^^-\'r'  +  (22) 

where  c!  =  dc/ds  may  be  non-zero  if  perturbations  are  present  (say,  if  c  depends  on  a  or  e).  In  the 
present  case  c  =  l/y^,  and  n  =  3/2,  so  (22)  can  be  solved  for  r", 

r"  =  —  ( \r^rr  +  .  (23) 


This  equation  involves  the  derivative  of  the  (scalar)  magnitude  r  which  is  easily  calculated, 

.  _  d^Jr  ■  r  _  r  ■  r 

dt  r  ’ 


so  that  the  second  derivative  may  be  rewritten, 


r  —  _  •  rjr  -h  r  r  j  .  [Zb) 

These  second  derivatives  are  then  integrated  using  the  Gauss- Jackson  and  summed  Adams  mid¬ 
corrector  formulas  to  find  the  position  and  velocity  at  each  of  the  8  points  surrounding  epoch.  The 
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integration  does  not  give  velocity  directly;  it  gives  r'  =  dr /ds.  So  a  conversion  must  be  made  to  find 
velocity, 


3 


(26) 


Before  the  force  model  can  be  re-evaluated,  the  time  at  each  point  must  be  found.  The  time  is  found 
be  integrating  a  seventh  differential  equation. 


t'  = 


1  3 


(27) 


using  the  summed  Adams  mid-corrector  formulas.  With  the  time  known  the  forces  are  evaluated  to 
compute  refined  estimates  of  the  accelerations,  these  accelerations  are  converted  into  s-space  second 
derivatives,  and  the  integration  is  performed  again  to  obtain  positions,  velocities,  and  times  at  the 
points.  This  process  repeats  until  the  accelerations  between  two  iterations  converges  to  a  prescribed 
tolerance.  The  initialization  procedure  for  s  integration  may  thus  be  summarized: 

Inititalization  for  s -integration 

1.  Convert  t  step  size  to  s  step  using  the  perigee  distance,  (18). 

2.  Use  /  and  g  series  to  calculate  8  positions  and  velocities  surrounding  epoch,  holding  the  epoch 
distance  constant  to  find  the  time,  (19). 

3.  Evaluate  9  accelerations  from  the  positions  and  velocities,  including  epoch. 

4.  Convert  the  accelerations  into  s  derivatives,  (25). 

5.  While  the  s  second  derivatives  have  not  converged: 

(a)  For  each  point  n  =  — 4 . . .  4,  n  7^  0: 

i.  Calculate  and  r'^,  using  Gauss- Jackson  and  summed  Adams  mid-corrector  formu¬ 
las. 

ii.  Convert  r'^  to  f*^,  (26). 

hi.  Calculate  the  time  at  point  n  by  integrating  (27)  with  the  summed  Adams  mid¬ 
corrector  formulas. 

iv.  Evaluate  the  acceleration  using  the  appropriate  force  model. 

V.  Convert  to  r",  (25). 

(b)  Test  convergence  of  r'f. 

After  the  integrator  has  completed  the  initialization  process,  it  goes  into  a  predictor-corrector 
cycle.  In  t-integration,  the  predictor  formulas  are  first  used  to  find  the  position  and  velocity  at 
the  next  point,  using  the  9  known  accelerations.  The  position  and  velocity  are  used  to  find  the 
acceleration  at  that  point,  and  that  acceleration  is  then  used  in  the  corrector  formulas  to  find  more 
accurate  values  of  the  position  and  velocity.  These  values  of  position  and  velocity  are  then  used  to 
find  a  more  accurate  value  of  the  acceleration  and  the  corrector  formulas  are  applied  again.  This 
process  repeats  until  the  position  and  velocity  converge  to  a  given  tolerance.  The  predictor  and 
corrector  cycles  for  t-integration  continue  from  steps  1-3  above  as  follows. 

Predict  (t -integration) 

4.  Calculate  and  using  Gauss- Jackson  and  summed  Adams  predictor  formulas. 

Evaluate  —  Correct  (t -integration) 

5.  Evaluate  the  acceleration 


6.  Increment  n. 
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7.  While  Vn  and  Vn  have  not  converged: 

(a)  Calculate  and  using  Gauss- Jackson  and  summed  Adams  corrector  formulas. 

(b)  Test  convergence  of  and  if  not  converged,  evaluate 

8.  Predict  next  time  step  (go  to  4). 

This  cycle  is  modified  for  s-integration.  When  the  independent  variable  is  s,  the  predictor 
(corrector)  does  not  directly  give  velocity,  it  gives  dr/ds,  so  the  velocity  must  be  found  using  (26). 
After  the  position  and  velocity  are  found  the  time  at  the  new  point  must  be  found  by  solving 
(27)  using  the  summed  Adams  predictor  (or  corrector)  formula.  After  the  forces  are  evaluated,  the 
accelerations  must  be  converted  to  second  s  derivatives  using  (25).  Thus  the  predictor-corrector 
modified  for  s-integration  continues  from  steps  1-5  above  as  follows: 

Predict  (s -integration) 

6.  Calculate  Vn+i  and  using  Gauss- Jackson  and  summed  Adams  predictor  formulas. 

7.  Convert  to  Vn+i,  (26). 

8.  Calculate  the  time  at  point  n  +  1  by  integrating  (27)  with  the  summed  Adams  predictor 
formula. 

Evaluate  —  Correct  (s -integration) 

10.  Evaluate  the  acceleration 

11.  Convert  Vn+i  to  (25). 

12.  Increment  n. 

13.  While  Vn  and  have  not  converged: 

(a)  Calculate  Vn  and  using  Gauss- Jackson  and  summed  Adams  corrector  formulas. 

(b)  Convert  r'^  to  (26). 

(c)  Calculate  the  time  at  point  n  by  integrating  (27)  with  the  summed  Adams  corrector 
formula. 

(d)  Test  convergence  of  and  if  not  converged,  evaluate  Vn,  and  convert  to  r". 

14.  Predict  next  time  step  (go  to  6). 

REGULARIZATION 

The  equations  of  motion  contain  the  distance  from  the  earth  center  to  the  satellite  in  the  denom¬ 
inator,  so  a  singularity  exists;  the  equations  are  not  regular.  The  equations  can  be  regularized  by 
introducing  two-body  conserved  quantities  (elements)  that  are  redundant,  such  as  the  energy  and  the 
Laplace  vector.  These  quantities  must  also  be  integrated  if  perturbations  are  present.  Schumacher 
has  suggested^  that  regularization  should  increase  the  accuracy  of  s-integration. 

The  equations  may  be  regularized  for  any  value  of  n  and  c  with  the  following  procedure  (Schu¬ 
macher,  Ref.  2,  pp.  17-24).  The  perturbed  equation  of  motion, 

r  +  ^r  =  P,  (28) 

^Private  communication. 
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where  P  is  the  perturbing  acceleration,  contains  a  singularity  at  r  =  0  because  of  the  reciprocal 
distance  term.  The  acceleration  r  in  (22)  can  be  replaced  using  the  equation  of  motion  (28), 


r"  —  nr  ^ r' r' +  =  c‘^r‘^'^P  +  c  ^c'r'.  (29) 

In  order  to  regularize  this  equation,  the  Keplerian  energy,  S,  and  the  Laplace  vector,  B,  must  be 
introduced, 

S=\r-r-f^,  (30) 

B  =  r  X  (r  X  r)  —  —r.  (31) 

r 

These  definitions  can  be  transformed  into  s-space  by  using  (21), 

£  =  ■  r')  - /ir-\  (32) 

B  =  X  {r  X  r')  —  fir~^r  (33) 

=  ■  r')  —  —  iir~^r.  (34) 

The  first  and  third  terms  of  the  last  equation  resemble  the  terms  in  the  energy  equation,  and  the 
second  term  resembles  the  non- regular  term  in  (29).  So  the  energy  and  the  Laplace  vector  can  be 
combined  to  regularize  (29), 


^{B  —  2Sr)  =  —nr  ^r'r' +  +  (n  —  l)c^/ir^^  (35) 

notice  that  the  term  has  been  deliberately  split  in  two.  Substituting  the  first  two  terms  on 

the  right  side  in  r"  equation  (29)  yields  a  regular  equation, 

r"  -  2Sc^nr^^~^r  -  {n  -  l)c^/ir^""“^r  =  -nc^r^'^~^B  +  c^r^'^P  +  c~^c'r',  (36) 

which,  because  of  the  third  term,  is  regular  for  n  =  1  or  n  >  3/2.  In  the  present  case  n  =  3/2  and 
c  =  so  (36)  becomes 

r"  =  —Err'  +  -r  —  ^rB  H — r^P,  (37) 

fi  2  2(1  11 

which  is  the  regularized  equivalent  of  (25).  When  perturbations  are  present  the  time  evolution  of 
energy  and  the  Laplace  vector  must  also  be  determined  by  integration.  To  find  the  derivative  of  the 
energy  first  take  a  time  derivative  of  (30), 

E  =  r  •  r  +  -^r.  (38) 

The  right  side  of  (38)  is  equivalent  to  taking  a  dot  product  of  r  with  the  equation  of  motion  (28), 


E  =  P  r.  (39) 

Both  sides  can  now  be  converted  to  s  derivatives, 

E'  =  P-r'.  (40) 

By  manipulating  the  equation  of  motion  (28),  the  time  derivative  of  the  Laplace  vector  (31)  can  be 
shown  to  be 


P  =  P  X  (r  X  r)  +  r{r  x  P) 

=  2r{P  ■  r)  -  r(P  ■  r)  -  P(r  ■  r). 


The  time  derivative  can  be  converted  to  an  s  derivative 


B'  =  2r{P  •  r')  -  r'{P  •  r)  -  P{r  •  r').  (43) 

Regularization  can  be  added  to  SpecialK  in  order  to  assess  the  effect  on  accuracy  and  computation 
speed.  After  the  forces  are  evaluated,  the  two  body  acceleration  is  subtracted  to  give  the  perturbation 
acceleration,  using  (28).  The  perturbation  acceleration  is  then  used  in  (37)  to  give  the  second 
derivative,  which  is  integrated  numerically  using  the  Gauss-Jackson  method.  The  energy  and  the 
Laplace  vector  are  also  found  by  numerically  integrating  (40)  and  (43),  respectively.  These  equations 
are  integrated  using  the  summed  Adams  method.  Equations  (30)  and  (31)  are  used  to  find  the  initial 
values  of  energy  and  the  Laplace  vector  in  the  initialization  routine.  The  procedure  can  be  described 
as  follows: 

Initialization 

1.  Convert  t  step  size  to  s  step  using  the  perigee  distance,  (18). 

2.  Find  the  initial  energy  and  Laplace  vector,  (30)  and  (31). 

3.  Use  /  and  g  series  to  calculate  8  positions  and  velocities  surrounding  epoch,  holding  the  epoch 
distance  constant  to  find  the  time,  (19). 

4.  Evaluate  9  accelerations  from  the  positions  and  velocities,  including  epoch. 

5.  Calculate  the  perturbing  acceleration  from  the  total  acceleration,  (28). 

6.  Calculate  the  s  derivatives,  (37). 

7.  While  the  s  second  derivatives  have  not  converged: 

(a)  For  each  point  n  =  —4 ...  4,  n  7^  0: 

i.  Calculate  and  r'^,  using  Gauss-Jackson  and  summed  Adams  mid-corrector  formu¬ 
las. 

ii.  Convert  r'^  to  (26). 

hi.  Find  the  energy  and  Laplace  vector  at  point  n  by  integrating  (40)  and  (43)  with  the 
summed  Adams  mid-corrector  formulas. 

iv.  Calculate  the  time  at  point  n  by  integrating  (27)  with  the  summed  Adams  mid¬ 
corrector  formulas. 

V.  Evaluate  the  acceleration  using  the  appropriate  force  model,  and  find  P. 
vi.  Calculate  r",  (37). 

(b)  Test  convergence  of  r". 

Predict 

8.  Calculate  Vn+i  and  using  Gauss-Jackson  and  summed  Adams  predictor  formulas. 

9.  Convert  to  Vn+i,  (26). 

10.  Calculate  the  energy,  Laplace  vector,  and  time  at  point  n  +  1  by  integrating  (40),  (43),  and 
(27)  with  the  summed  Adams  predictor  formula. 

Evaluate  —  Correct 

11.  Evaluate  the  acceleration  r^+i,  and  find  P. 

12.  Calculate  (37). 

13.  Increment  n. 
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14.  While  Vn  and  Vn  have  not  converged: 

(a)  Calculate  and  using  Gauss- Jackson  and  summed  Adams  corrector  formulas. 

(b)  Convert  r'^  to  rn,  (26). 

(c)  Calculate  the  energy,  Laplace  vector,  and  time  at  point  n  by  integrating  (40),  (43),  and 
(27)  with  the  summed  Adams  corrector  formula. 

(d)  Test  convergence  of  and  if  not  converged,  evaluate  P,  and  calculate  r". 

15.  Predict  next  time  step  (go  to  8). 

COMPARISONS  WITH  t-INTEGRATION 

A  proposal  to  replace  an  existing  t-integrator  with  an  s-integrator  for  some  or  all  orbits  should 
weigh  two  considerations:  the  relative  computation  time  and  the  relative  accuracy.  Two  tests  of 
t-integration  and  s-integration  are  considered  in  order  to  assess  these  factors: 

•  For  testing  without  perturbations,  the  reference  values  are  taken  from  the  analytic  two-body 
solutions. 

•  For  testing  with  perturbations  (24  x  24  WGS-84  geopotential,  Jacchia  70  drag  model,  and 
lunar  and  solar  perturbations),  the  reference  values  are  taken  by  taking  the  final  point  of  the 
integration  and  integrating  backwards  to  epoch.  This  forward-backward  test  gives  a  rough 
indication  of  the  integration  error. 

In  all  tests,  a  metric  for  integration  accuracy  is  defined  using  an  error  ratio  defined  in  terms  of 
the  RMS  error  of  the  integration  (Ref.  13).  First  define  position  and  velocity  errors  as 


Ar  I  ^computed  ^reference  |? 
An  =  lUcomputed  '^reference  |- 


The  RMS  position  error  can  be  calculated. 


A^rms  = 


(44) 

(45) 


(46) 


with  a  similar  equation  for  RMS  velocity  error.  The  RMS  position  error  is  normalized  by  the  apogee 
distance  and  the  number  of  orbits  to  find  the  position  error  ratio, 

_  ArRMS 
ACrbits 

The  velocity  error  ratio  is  found  by  normalizing  by  the  perigee  speed  and  the  number  or  orbits, 

_  A-CrmS 
ACrbits 

Six  orbits  are  considered,  combinations  of  eccentricities  of  0.0,  0.25,  and  0.75  and  perigee  heights 
of  300  km,  and  1000  km.  All  of  the  orbits  have  an  inclination  of  40°,  and  a  ballistic  coefficient  of 
0.01  The  epoch  is  2001-10-01  00:00:00  UTC.  Error  ratios  for  position  and  velocity  are  found 

for  each  orbit  with  and  without  perturbations,  using  t-integration,  s-integration,  and  s-integration 
with  regularization.  The  error  ratio  for  s-integration  with  perturbations  using  t-integration  as  the 
reference  orbit  is  found  with  and  without  regularization.  Each  integration  is  for  three  days  (72 
hours),  with  a  30  second  time  step  for  t-integration,  and  the  corresponding  s  step  given  by  (18)  for 
s-integration.  The  ephemeris  is  generated  at  one  minute  time  intervals. 


(48) 
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In  cases  where  the  step  size  used  does  not  match  the  desired  output  points,  SpecialK  uses  an 
interpolator.  For  these  tests,  the  t-integration  output  points  directly  match  the  integration  points, 
so  no  interpolation  is  needed.  Interpolation  is  still  needed,  however,  for  s-integration,  where  the 
integration  time  does  not  occur  at  the  step  points.  The  interpolator  normally  used  in  SpecialK  is  a 
fifth-order  interpolator  for  position,  and  a  fourth-order  interpolator  for  velocity;  however,  in  order 
to  get  a  fair  comparison  of  t-  and  s-integration,  an  eighth-order  Lagrange  interpolator  developed 
for  these  tests  is  used  for  the  s-integration  runs.  This  is  the  same  order  as  the  eighth  order  Gauss- 
Jackson  integrator,  and  should  insure  no  loss  of  accuracy  during  the  s-integration  runs. 

Table  1  shows  the  results  of  the  two-body  test.  Ephemeris  is  generated  by  SpecialK  with  all 
the  perturbations  turned  off,  and  compared  to  the  analytic  solution.  The  ephemeris  generated  by 
t-integration  is  more  accurate,  most  notably  for  the  circular  orbit.  As  the  eccentricity  is  increased, 
the  accuracy  of  s-integration  approaches  that  of  t-integration.  Regularizing  the  s-integration  does 
not  show  a  significant  improvement  in  accuracy. 


Table  1:  Error  ratios  for  the  two-body  test  relative  to  an  analytic  solution  for  t-integration, 
integration,  and  s-integration  with  regularization. 


Test  Gase 

Error 

Ratio 

e 

hp  (km) 

t 

s 

s-reg 

0 

300 

pos 

8.40 

X 

10- 

-ir 

8.94 

X 

10- 

-12 

8.94 

X 

10- 

-12 

vel 

8.40 

X 

10- 

-17 

8.94 

X 

10- 

-12 

8.94 

X 

10- 

-12 

0 

1000 

pos 

7.36 

X 

10- 

-17 

4.33 

X 

10- 

-11 

4.33 

X 

10- 

-11 

vel 

7.36 

X 

10- 

-17 

4.33 

X 

10- 

-11 

4.33 

X 

10- 

-11 

0.25 

300 

pos 

8.05 

X 

10- 

-16 

1.47 

X 

10- 

-13 

1.40 

X 

10- 

-13 

vel 

8.60 

X 

10- 

-16 

1.57 

X 

10- 

-13 

1.49 

X 

10- 

-13 

0.25 

1000 

pos 

8.25 

X 

10- 

-17 

1.43 

X 

10- 

-13 

1.45 

X 

10- 

-13 

vel 

8.75 

X 

10- 

-17 

1.52 

X 

10- 

-13 

1.54 

X 

10- 

-13 

0.75 

300 

pos 

6.76 

X 

10- 

-15 

1.55 

X 

10- 

-14 

1.95 

X 

10- 

-14 

vel 

1.25 

X 

10- 

-14 

3.17 

X 

10- 

-14 

3.73 

X 

10- 

-14 

0.75 

1000 

pos 

1.04 

X 

10- 

-15 

1.14 

X 

10- 

-13 

1.07 

X 

10- 

-13 

vel 

2.27 

X 

10- 

-15 

2.42 

X 

10- 

-13 

2.27 

X 

10- 

-13 

s- 


Table  2  shows  the  results  of  testing  with  perturbations.  Each  orbit  is  propagated  forward  3  days 
with  a  24  X  24  WGS-84  geopotential,  Jacchia  70  drag  model,  and  lunar  and  solar  perturbations.  The 
final  value  of  the  propagation  is  then  used  to  propagate  backwards  to  the  original  epoch.  The  error 
ratio  shown  is  the  difference  between  the  forward  and  backward  propagations.  Again  s-integration 
has  a  higher  error  ratio  than  t-integration,  with  a  larger  difference  for  the  1000  km  orbits.  In  the 
presence  of  higher  drag,  t-integration  and  s-integration  both  lose  accuracy  and  have  closer  error 
ratios.  Again,  s-integration  with  regularization  has  similar  error  ratios  to  s-integration  without 
regularization,  with  some  improvement  showing  at  the  higher  eccentricity. 

Table  3  shows  the  difference  between  the  t-integration  and  s-integration  forward  propagations, 
as  well  as  the  difference  between  the  t-integration  and  s-integration  with  regularization  forward 
propagations.  The  higher  perigee  orbits  have  a  closer  agreement  between  t-  and  s-integration,  which 
is  consistent  with  the  results  from  Table  2.  Again  the  results  for  s-integration  with  and  without 
regularization  are  similar. 

A  speed  test  is  also  performed  for  each  orbit.  The  amount  of  time,  in  seconds,  to  propagate  30 
days  from  epoch  is  shown  for  each  integration  type  in  Table  4.  The  test  is  performed  on  an  SGI 
Origin  200,  and  the  time  shown  is  the  user  time.  The  300  km  circular  orbit  decays  after  21  days,  so 
the  time  shown  is  for  a  20  day  run.  Eor  the  circular  orbits,  the  t-  and  s-integration  take  the  same 
amount  of  time,  because  there  is  no  savings  in  integration  steps  for  s-integration  for  a  circular  orbit. 
Eor  the  0.25  eccentricity  orbits  there  is  a  1.5:1  time  savings,  and  for  the  0.75  eccentricity  orbits  there 
is  a  6:1  savings.  The  s-integration  with  regularization  has  only  a  marginal  cost  in  time  over  the 
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Table  2:  Error  ratios  for  test  with  perturbations,  backward  integration  relative  to  forward. 


Test  Case 

Error 

Ratio 

e 

hp  (km) 

t 

s 

s-reg 

0 

300 

pos 

4.93 

X 

10- 

-9 

1.69 

X 

10- 

-8 

1.69 

X 

10- 

-8 

vel 

4.93 

X 

10- 

-9 

1.69 

X 

10- 

-8 

1.69 

X 

10- 

-8 

0 

1000 

pos 

3.38 

X 

10- 

-12 

1.18 

X 

10- 

-9 

1.18 

X 

10- 

-9 

vel 

3.38 

X 

10- 

-12 

1.19 

X 

10- 

-9 

1.18 

X 

10- 

-9 

0.25 

300 

pos 

4.17 

X 

10- 

-10 

9.38 

X 

10- 

-9 

9.31 

X 

10- 

-9 

vel 

4.44 

X 

10- 

-10 

1.00 

X 

10- 

-8 

9.93 

X 

10- 

-9 

0.25 

1000 

pos 

2.17 

X 

10- 

-13 

2.72 

X 

10- 

-10 

2.80 

X 

10- 

-10 

vel 

2.32 

X 

10- 

-13 

2.90 

X 

10- 

-10 

2.99 

X 

10- 

-10 

0.75 

300 

pos 

5.78 

X 

10- 

-9 

1.61 

X 

10- 

-8 

1.56 

X 

10- 

-8 

vel 

1.23 

X 

10- 

-8 

3.38 

X 

10- 

-8 

3.29 

X 

10- 

-8 

0.75 

1000 

pos 

3.61 

X 

10- 

-12 

1.92 

X 

10- 

-10 

1.18 

X 

10- 

-10 

vel 

7.62 

X 

10- 

-12 

4.07 

X 

10- 

-10 

2.48 

X 

10- 

-10 

Table  3:  Error  ratio  of  s-integration  and  s-integration  using  regularization  with  t-integration  as  the 
reference,  with  perturbations. 


Test  Case 
e  hp  (km) 

Position  Error  Ratio 

s  s-reg 

Velocity  Error  Ratio 
s  s-reg 

0 

300 

1.24  X  10-« 

1.25  X  10-^ 

1.24  X  10-« 

1.25  X  10-« 

0 

1000 

4.00  X  10-^° 

4.00  X  10-^9 

4.00  X  10-^9 

4.00  X  10-19 

0.25 

300 

1.09  X  10-9 

1.14  X  10-9 

1.17  X  10-9 

1.22  X  10-9 

0.25 

1000 

9.19  X  10-^1 

9.40  X  10-1^ 

9.74  X  10-^1 

9.97  X  10-11 

0.75 

300 

4.42  X  10-s 

4.30  X  10-s 

9.41  X  10-s 

9.16  X  10-s 

0.75 

1000 

9.12  X  10-^1 

4.46  X  10-1^ 

1.95  X  10-^9 

9.22  X  10-11 

16 


non-regularized  s-integration. 


Table  4:  User  time  on  an  SGI  Origin  200  to  integrate  30  {20)  days  with  perturbations. 


Test  Gase 
e  hp  (km) 

Time  for  30  Day  Run  (sec) 
t  s  s-reg 

0 

300 

21 

21 

0 

1000 

31 

31 

32 

0.25 

300 

29 

20 

21 

0.25 

1000 

29 

20 

20 

0.75 

300 

28 

4.7 

4.8 

0.75 

1000 

28 

4.6 

4.7 

The  previous  results  show  that  when  the  step  size  for  s-integration  is  chosen  so  that  at  perigee 
it  uses  the  same  step  size  as  t-integration,  s-integration  is  faster  for  integrating  elliptical  orbits  but 
is  not  as  accurate.  This  conclusion  leads  to  another  study,  where  the  step  size  is  chosen  to  give 
comparable  accuracy  between  t-integration  and  s-integration.  In  this  set  of  tests,  the  step  size  for 
each  test  case  is  found  that  gives  an  error  ratio  of  approximately  1  x  10“^  in  the  forward-backward 
test.  This  step  size  is  then  used  in  a  speed  test,  once  again  propagating  30  days  forward  from  epoch 
with  perturbations  turned  on.  This  test  is  performed  to  compare  t-integration  to  s-integration 
without  regularization;  the  previous  results  indicate  that  regularization  does  not  yield  a  signihcant 
gain. 

Table  5  shows  the  results  for  orbits  with  a  1000  km  perigee  height.  Again,  all  of  the  test  cases 
have  an  inclination  of  40°,  and  a  ballistic  coefficient  of  0.01  The  table  shows  the  step  size  that 

gives  an  error  ratio  of  1  x  lO""^  in  the  forward-backward  test,  and  the  time  to  propagate  forward 
30  days  using  that  step  size,  for  both  t-integration  and  s-integration.  The  table  also  gives  the 
ratio  of  the  run-time  for  t-integration  to  the  run-time  for  s-integration.  When  the  ratio  is  greater 
than  1,  s-integration  is  preferable.  Table  5  shows  that  s-integration  is  preferable  for  e  >  0.15.  At 
eccentricities  of  0.20  and  0.25,  the  step  size  used  in  s-integration  is  actually  larger  than  the  step 
size  for  t-integration,  which  indicates  that  for  these  step  sizes  s-integration  is  more  accurate  than 
t-integration.  This  is  contrary  to  the  results  shown  in  Table  2,  and  this  may  mean  that  at  the  smaller 
step  sizes  used  in  Table  2,  s-integration  is  affected  by  round-off  error. 


Table  5:  User  time  on  an  SGI  Origin  200  to  integrate  30  days  with  perturbations,  with  the  step  size 
shown.  Step  size  chosen  to  give  an  error  ratio  of  1  x  10“^.  Perigee  height  is  1000  km. 


e 

Step  Size 
t  s 

Time  for  30  Day  Run  (sec) 
t  s 

Ratio 

0 

50 

30 

18.3 

31.0 

0.59 

0.05 

55 

54 

16.7 

16.4 

1.0 

0.10 

64 

58 

14.2 

14.0 

1.0 

0.15 

73 

72 

12.2 

10.3 

1.2 

0.20 

70 

74 

12.6 

9.17 

1.4 

0.25 

68 

80 

12.8 

7.80 

1.6 

0.50 

61 

61 

14.0 

5.80 

2.4 

0.75 

65 

51 

13.1 

2.94 

4.5 

The  results  for  a  perigee  height  of  300  km  are  shown  in  Table  6.  In  this  case,  s-integration  is 
preferable  for  e  >  0.30.  The  primary  reason  why  these  results  differ  from  those  of  1000  km  is  that  the 
difference  between  the  t-integration  and  s-integration  step  sizes  needed  to  maintain  a  1  x  10“^  error 
ratio  is  more  significant.  In  the  presence  of  higher  drag,  both  t-integration  and  s-integration  require 
a  smaller  step  size  than  needed  for  the  1000  km  case.  However  the  step  size  has  to  be  reduced  more 
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dramatically  for  s-integration.  Again  this  indicates  that  s-integration  is  being  affected  by  round-off 
error  more  severely  than  t-integration.  In  fact,  while  searching  for  the  step  sizes  needed  for  this 
test,  we  found  that  in  some  cases  lowering  the  step  size  increased  the  error  ratio.  That  may  also 
be  due  to  the  nature  of  the  forward-backward  test,  however.  If  s-integration  is  being  affected  by 
round-off  error  more  severely  than  t-integration,  the  error  is  most  likely  coming  from  integrating  the 
time  equation,  which  is  only  needed  in  s-integration.  More  study  is  needed  to  determine  if  the  error 
from  integrating  the  time  equation  can  be  reduced. 


Table  6:  User  time  on  an  SGI  Origin  200  to  integrate  30  {20)  days  from  epoch  with  perturbations, 
with  the  step  size  shown.  Step  size  chosen  to  give  an  error  ratio  of  1  x  10“ Perigee  height  is  300 
km.  _ 


e 

Step  Size 
t  s 

Time  for  30  Day  Run  (sec) 
t  s 

Ratio 

0 

6 

4 

99.6 

153 

0.65 

0.15 

50 

26 

18.1 

27.9 

0.65 

0.25 

32 

15 

27.0 

38.9 

0.69 

0.30 

17 

16 

50.0 

33.0 

1.5 

0.35 

32 

30 

26.7 

16.2 

1.6 

0.40 

50 

36 

17.2 

12.2 

1.4 

0.50 

50 

36 

17.1 

9.55 

1.8 

0.75 

29 

10 

28.8 

13.0 

2.2 

CONCLUSION  AND  FURTHER  CONSIDERATION 

It  is  possible  to  accelerate  the  numerical  integration  of  orbits  without  loss  of  accuracy  by  transform¬ 
ing  the  independent  variable  from  time  (mean  anomaly)  to  another  variable  s  by  the  generalized 
Sundman  transformation  dt  =  r^ds  for  some  exponent  n.  Existing  literature  indicates  n  =  3/2  is 
optimal  for  this  purpose.  For  circular  orbits,  the  two  integration  methods  s  and  t  have  exactly  the 
same  number  of  points;  for  eccentric  orbits,  s-integration  has  fewer  points  on  an  orbit  for  equiv¬ 
alent  accuracy.  One  might  be  tempted  to  conclude  that  s-integration  is  always  indicated  because 
it  is  never  worse  than  t-integration.  This  is  not  so,  however.  A  practical  s-integration  requires 
the  integration  of  a  seventh  differential  equation,  the  time  equation,  and  also  requires  additional 
transformations  to  and  from  time.  These  have  the  effect  of  costing  a  small  but  measurable  amount 
of  computation  time.  Furthermore,  the  introduction  of  perturbations  causes  proportionately  more 
inaccuracies  than  two-body  motion  because  of  the  less-frequent  force  evaluation.  For  sufficiently  low 
eccentricities,  t-integration  is  faster  for  a  given  accuracy.  A  two-body  analysis  alone  is  not  sufficient 
to  determine,  for  a  given  orbit,  which  method  is  faster. 

We  sought  to  determine  a  threshold  eccentricity  above  which  s-integration  is  to  be  preferred 
because  of  decreased  computation  time  for  equivalent  accuracy.  It  turns  out  this  threshold  is  de¬ 
pendent  on  other  orbit  parameters,  specifically,  perigee  height.  Based  on  a  study  of  the  s  and  t 
integrators  in  the  SpecialK  code  with  step  size  set  so  that  errors  are  comparable  and  a  realistic  set 
of  perturbation  forces,  we  have  found  that  for  a  perigee  height  of  1000  km,  the  threshold  is  e  0.15. 
At  a  lower  perigee  of  300  km,  the  threshold  is  e  ~  0.30.  We  believe  the  amount  of  drag  the  satellite 
experiences  affects  the  round-off  error  and  contributes  to  the  determination  of  this  threshold. 

We  found  that  regularization  of  the  differential  equations  results  in  a  minor  improvement  in 
accuracy  for  eccentricities  higher  than  the  threshold  and  does  not  justify  the  implementation  or 
computation  time. 

Nacozy  (Ref.  17)  studied  the  time  transformation  equation  (2)  with  the  intent  of  finding  alternate 
forms  that  may  improve  accuracy  of  the  numerical  integration  of  that  equation.  He  suggested  several 
reformulations  of  the  time  transformation  equation  that  results  in  accuracy  improvements  of  the 
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numerical  integration  of  the  time  equation,  including  one  that  appears  to  be  optimal  for  n=3/2.  His 
test  included  only  the  two-body  force  and  the  perturbation  on  a  satellite  with  e  =  0.3,  and  no 
computation  time  cost  is  assessed.  However,  this  work  may  be  a  fruitful  area  for  further  research  in 
improving  the  s-integration  results. 

Another  accuracy-related  consideration  that  we  have  not  addressed  is  stability  in  the  Lyupanov 
sense:  how  much  does  a  given  small  error  in  the  orbital  parameters  contribute  to  error  in  the  central 
angle,  and  thus  to  in-track  position  error,  over  a  given  period  of  time?  Baumgarte  (Ref.  18)  discussed 
this  problem  and  determined  that  the  KS  transformation  (n  =  1)  improves  orbit  integration  in  this 
regard.  This  problem  can  be  important  because,  while  one  integrator  may  have  the  same  or  even 
greater  error  ratio  than  another,  it  may  still  be  more  stable  in  the  Lyapunov  sense,  and  would  keep 
the  in- track  errors  more  confined.  Keeping  the  in- track  errors  confined  increases  the  time  span  an 
orbit  propagation  remains  accurate. 
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